module physpkg

  !-----------------------------------------------------------------------
  ! Purpose:
  !
  ! Provides the interface to CAM physics package
  !
  ! Revision history:
  ! Aug  2005,  E. B. Kluzek,  Creation of module from physpkg subroutine
  ! 2005-10-17  B. Eaton       Add contents of inti.F90 to phys_init().  Add
  !                            initialization of grid info in phys_state.
  ! Nov 2010    A. Gettelman   Put micro/macro physics into separate routines
  !-----------------------------------------------------------------------

  use shr_kind_mod,     only: r8 => SHR_KIND_R8
  use spmd_utils,       only: masterproc, iam
  use physconst,        only: latvap, latice, rh2o
  use physics_types,    only: physics_state, physics_tend, physics_state_set_grid, &
                              physics_ptend, physics_tend_init, physics_update,    &
                              physics_type_alloc, physics_ptend_dealloc,           &
                              physics_state_alloc, physics_state_dealloc,          &
                              physics_tend_alloc, physics_tend_dealloc
  use phys_grid,        only: get_ncols_p
  use phys_gmean,       only: gmean_mass
  use ppgrid,           only: begchunk, endchunk, pcols, pver, pverp, psubcols
  use constituents,     only: pcnst, cnst_name, cnst_get_ind
  use camsrfexch,       only: cam_out_t, cam_in_t
  use cam_control_mod,  only: ideal_phys, adiabatic
  use phys_control,     only: phys_do_flux_avg, phys_getopts, waccmx_is
  use flux_avg,         only: flux_avg_init
  use cam_logfile,      only: iulog
  use camsrfexch,       only: cam_export

  use modal_aero_calcsize,    only: modal_aero_calcsize_init, modal_aero_calcsize_diag, modal_aero_calcsize_reg
  use modal_aero_wateruptake, only: modal_aero_wateruptake_init, modal_aero_wateruptake_dr, modal_aero_wateruptake_reg

  implicit none

  private

  save

  ! Public methods
  public phys_register ! was initindx  - register physics methods
  public phys_init     ! Public initialization method
  public phys_run1     ! First phase of the public run method
  public phys_run2     ! Second phase of the public run method
  public phys_final    ! Public finalization method

  ! Physics package options
  character(16) shallow_scheme
  character(16) macrop_scheme
  character(16) microp_scheme
  character(16) subcol_scheme
  integer       cld_macmic_num_steps      ! Number of macro/micro substeps
  logical       do_clubb_sgs
  logical       use_subcol_microp         ! if true, use subcolumns in microphysics
  logical       state_debug_checks        ! Debug physics_state.
  logical       clim_modal_aero           ! climate controled by prognostic or prescribed modal aerosols
  logical       prog_modal_aero           ! Prognostic modal aerosols present

  !  Physics buffer index
  integer :: teout_idx          = 0
  integer :: landm_idx          = 0
  integer :: sgh_idx            = 0
  integer :: sgh30_idx          = 0
  integer :: qini_idx           = 0
  integer :: cldliqini_idx      = 0
  integer :: cldiceini_idx      = 0
  integer :: prec_str_idx       = 0
  integer :: snow_str_idx       = 0
  integer :: prec_sed_idx       = 0
  integer :: snow_sed_idx       = 0
  integer :: prec_pcw_idx       = 0
  integer :: snow_pcw_idx       = 0
  integer :: prec_dp_idx        = 0
  integer :: snow_dp_idx        = 0
  integer :: prec_sh_idx        = 0
  integer :: snow_sh_idx        = 0
  integer :: dlfzm_idx          = 0     ! detrained convective cloud water mixing ratio.

contains

  subroutine phys_register()

    use cam_abortutils      , only: endrun
    use physics_buffer      , only: pbuf_init_time
    use physics_buffer      , only: pbuf_add_field, dtype_r8, pbuf_register_subcol
    use shr_kind_mod        , only: r8 => shr_kind_r8
    use spmd_utils          , only: masterproc
    use constituents        , only: pcnst, cnst_add, cnst_chk_dim, cnst_name
    use cam_control_mod     , only: moist_physics
    use chemistry           , only: chem_register
    use cloud_fraction      , only: cldfrc_register
    use rk_stratiform       , only: rk_stratiform_register
    use microp_driver       , only: microp_driver_register
    use microp_aero         , only: microp_aero_register
    use macrop_driver       , only: macrop_driver_register
    use clubb_intr          , only: clubb_register_cam
    use conv_water          , only: conv_water_register
    use physconst           , only: mwdry, cpair, mwh2o, cpwv
    use check_energy        , only: check_energy_register
    use carma_intr          , only: carma_register
    use cam3_aero_data      , only: cam3_aero_data_on, cam3_aero_data_register
    use cam3_ozone_data     , only: cam3_ozone_data_on, cam3_ozone_data_register
    use ghg_data            , only: ghg_data_register
    use vertical_diffusion  , only: vd_register
    use convect_deep        , only: convect_deep_register
    use convect_shallow     , only: convect_shallow_register
    use radiation           , only: radiation_register
    use co2_cycle           , only: co2_register
    use flux_avg            , only: flux_avg_register
    use iondrag             , only: iondrag_register
    use waccmx_phys_intr    , only: waccmx_phys_ion_elec_temp_reg
    use string_utils        , only: to_lower
    use prescribed_ozone    , only: prescribed_ozone_register
    use prescribed_volcaero , only: prescribed_volcaero_register
    use prescribed_strataero, only: prescribed_strataero_register
    use prescribed_aero     , only: prescribed_aero_register
    use prescribed_ghg      , only: prescribed_ghg_register
    use sslt_rebin          , only: sslt_rebin_register
    use aoa_tracers         , only: aoa_tracers_register
    use aircraft_emit       , only: aircraft_emit_register
    use cam_diagnostics     , only: diag_register
    use cloud_diagnostics   , only: cloud_diagnostics_register
    use cospsimulator_intr  , only: cospsimulator_intr_register
    use rad_constituents    , only: rad_cnst_get_info ! Added to query if it is a modal aero sim or not
    use subcol              , only: subcol_register
    use subcol_utils        , only: is_subcol_on, subcol_get_scheme
    use dyn_comp            , only: dyn_register

    integer m, nmodes

    ! Get physics options
    call phys_getopts(shallow_scheme_out          =shallow_scheme      , &
                      macrop_scheme_out           =macrop_scheme       , &
                      microp_scheme_out           =microp_scheme       , &
                      cld_macmic_num_steps_out    =cld_macmic_num_steps, &
                      do_clubb_sgs_out            =do_clubb_sgs        , &
                      use_subcol_microp_out       =use_subcol_microp   , &
                      state_debug_checks_out      =state_debug_checks  )

    subcol_scheme = subcol_get_scheme()

    ! Initialize dyn_time_lvls
    call pbuf_init_time()

    ! Register the subcol scheme
    call subcol_register()

    ! Register water vapor.
    ! ***** N.B. ***** This must be the first call to cnst_add so that
    !                  water vapor is constituent 1.
    if (moist_physics) then
      call cnst_add('Q', mwh2o, cpwv, 1.0e-12_r8, m, &
        longname='Specific humidity', readiv=.true., is_convtran1=.true.)
    else
      call cnst_add('Q', mwh2o, cpwv, 0.0_r8, m, &
        longname='Specific humidity', readiv=.false., is_convtran1=.true.)
    end if

    ! Topography file fields.
    call pbuf_add_field('LANDM',     'global',  dtype_r8, [pcols],      landm_idx)
    call pbuf_add_field('SGH',       'global',  dtype_r8, [pcols],      sgh_idx)
    call pbuf_add_field('SGH30',     'global',  dtype_r8, [pcols],      sgh30_idx)

    ! Fields for physics package diagnostics
    call pbuf_add_field('QINI',      'physpkg', dtype_r8, [pcols,pver], qini_idx)
    call pbuf_add_field('CLDLIQINI', 'physpkg', dtype_r8, [pcols,pver], cldliqini_idx)
    call pbuf_add_field('CLDICEINI', 'physpkg', dtype_r8, [pcols,pver], cldiceini_idx)

    ! check energy package
    call check_energy_register()

    ! If using a simple physics option (e.g., held_suarez, adiabatic),
    ! the normal CAM physics parameterizations are not called.
    if (moist_physics) then
      ! register fluxes for saving across time
      if (phys_do_flux_avg()) call flux_avg_register()

      call cldfrc_register()

      ! cloud water
      if (microp_scheme == 'RK') then
        call rk_stratiform_register()
      else if (microp_scheme == 'MG') then
        if (.not. do_clubb_sgs) call macrop_driver_register()
        call microp_aero_register()
        call microp_driver_register()
      end if

      ! Register CLUBB_SGS here
      if (do_clubb_sgs) call clubb_register_cam()

      call pbuf_add_field('PREC_STR', 'physpkg', dtype_r8, [pcols], prec_str_idx)
      call pbuf_add_field('SNOW_STR', 'physpkg', dtype_r8, [pcols], snow_str_idx)
      call pbuf_add_field('PREC_PCW', 'physpkg', dtype_r8, [pcols], prec_pcw_idx)
      call pbuf_add_field('SNOW_PCW', 'physpkg', dtype_r8, [pcols], snow_pcw_idx)
      call pbuf_add_field('PREC_SED', 'physpkg', dtype_r8, [pcols], prec_sed_idx)
      call pbuf_add_field('SNOW_SED', 'physpkg', dtype_r8, [pcols], snow_sed_idx)
      if (is_subcol_on()) then
        call pbuf_register_subcol('PREC_STR', 'phys_register', prec_str_idx)
        call pbuf_register_subcol('SNOW_STR', 'phys_register', snow_str_idx)
        call pbuf_register_subcol('PREC_PCW', 'phys_register', prec_pcw_idx)
        call pbuf_register_subcol('SNOW_PCW', 'phys_register', snow_pcw_idx)
        call pbuf_register_subcol('PREC_SED', 'phys_register', prec_sed_idx)
        call pbuf_register_subcol('SNOW_SED', 'phys_register', snow_sed_idx)
      end if

      ! Who should add FRACIS?
      ! -- It does not seem that aero_intr should add it since FRACIS is used in convection
      !     even if there are no prognostic aerosols ... so do it here for now
      call pbuf_add_field('FRACIS', 'physpkg', dtype_r8, [pcols,pver,pcnst], m)

      call conv_water_register()

      ! Determine whether its a 'modal' aerosol simulation or not
      call rad_cnst_get_info(0, nmodes=nmodes)
      clim_modal_aero = (nmodes > 0)

      if (clim_modal_aero) then
        call modal_aero_calcsize_reg()
        call modal_aero_wateruptake_reg()
      end if

      ! register chemical constituents including aerosols ...
      call chem_register()

      ! co2 constituents
      call co2_register()

      ! register data model ozone with pbuf
      if (cam3_ozone_data_on) then
        call cam3_ozone_data_register()
      end if
      call prescribed_volcaero_register()
      call prescribed_strataero_register()
      call prescribed_ozone_register()
      call prescribed_aero_register()
      call prescribed_ghg_register()
      call sslt_rebin_register

      ! CAM3 prescribed aerosols
      if (cam3_aero_data_on) then
        call cam3_aero_data_register()
      end if

      ! register various data model gasses with pbuf
      call ghg_data_register()

      ! carma microphysics
      call carma_register()

      ! Register iondrag variables with pbuf
      call iondrag_register()

      ! Register ionosphere variables with pbuf if mode set to ionosphere
      if (waccmx_is('ionosphere')) then
        call waccmx_phys_ion_elec_temp_reg()
      end if

      call aircraft_emit_register()

      ! deep convection
      call convect_deep_register()

      ! shallow convection
      call convect_shallow_register()

      ! radiation
      call radiation_register()
      call cloud_diagnostics_register()

      ! COSP
      call cospsimulator_intr_register()

      ! vertical diffusion
      call vd_register()
    else
      ! held_suarez/adiabatic physics option should be in simple_physics
      call endrun('phys_register: moist_physics configuration error')
    end if

    ! Register diagnostics PBUF
    call diag_register()

    ! Register age of air tracers
    call aoa_tracers_register()

    call dyn_register()

    ! All tracers registered, check that the dimensions are correct
    call cnst_chk_dim()

    ! ***NOTE*** No registering constituents after the call to cnst_chk_dim.

  end subroutine phys_register

  subroutine phys_inidat(cam_out, pbuf2d)

    use cam_abortutils,      only: endrun
    use physics_buffer,      only: pbuf_get_index, pbuf_get_field, physics_buffer_desc, pbuf_set_field, dyn_time_lvls
    use cam_initfiles,       only: initial_file_get_id, topo_file_get_id
    use cam_grid_support,    only: cam_grid_check, cam_grid_id
    use cam_grid_support,    only: cam_grid_get_dim_names
    use pio,                 only: file_desc_t
    use ncdio_atm,           only: infld
    use dycore,              only: dycore_is
    use polar_avg,           only: polar_average
    use short_lived_species, only: initialize_short_lived_species
    use cam_control_mod,     only: aqua_planet
    use waccmx_phys_intr,    only: waccmx_phys_ion_elec_temp_inidat

    type(cam_out_t), intent(inout) :: cam_out(begchunk:endchunk)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)

    integer lchnk, m, n, i, k, ncol
    type(file_desc_t), pointer :: fh_ini, fh_topo
    character(8) fieldname
    real(r8), pointer :: tptr(:,:), tptr_2(:,:), tptr3d(:,:,:), tptr3d_2(:,:,:)
    real(r8), pointer :: qpert(:,:)

    character(11), parameter :: sub = 'phys_inidat'
    integer tpert_idx, qpert_idx, pblh_idx

    logical :: found, found2=.false.
    integer ierr
    character(8) dim1name, dim2name
    integer ixcldice, ixcldliq
    integer grid_id  ! grid ID for data mapping

    nullify(tptr, tptr_2, tptr3d, tptr3d_2)

    fh_ini  => initial_file_get_id()
    fh_topo => topo_file_get_id()

    ! dynamics variables are handled in dyn_init - here we read variables needed for physics
    ! but not dynamics

    grid_id = cam_grid_id('physgrid')
    if (.not. cam_grid_check(grid_id)) then
      call endrun(trim(sub)//': Internal error, no "physgrid" grid')
    end if
    call cam_grid_get_dim_names(grid_id, dim1name, dim2name)

    allocate(tptr(1:pcols,begchunk:endchunk))

    if (associated(fh_topo) .and. .not. aqua_planet) then
      call infld('SGH', fh_topo, dim1name, dim2name, 1, pcols, begchunk, endchunk, tptr, found, gridname='physgrid')
      if (.not. found) call endrun('ERROR: SGH not found on topo file')

      call pbuf_set_field(pbuf2d, sgh_idx, tptr)

      allocate(tptr_2(1:pcols,begchunk:endchunk))
      call infld('SGH30', fh_topo, dim1name, dim2name, 1, pcols, begchunk, endchunk, tptr_2, found, gridname='physgrid')
      if (found) then
        call pbuf_set_field(pbuf2d, sgh30_idx, tptr_2)
      else
        if (masterproc) write(iulog,*) 'Warning: Error reading SGH30 from topo file.'
        if (masterproc) write(iulog,*) 'The field SGH30 will be filled using data from SGH.'
        call pbuf_set_field(pbuf2d, sgh30_idx, tptr)
      end if

      deallocate(tptr_2)

      call infld('LANDM_COSLAT', fh_topo, dim1name, dim2name, 1, pcols, begchunk, endchunk, tptr, found, gridname='physgrid')
      if (.not.found) call endrun(' ERROR: LANDM_COSLAT not found on topo dataset.')

      call pbuf_set_field(pbuf2d, landm_idx, tptr)
    else
      call pbuf_set_field(pbuf2d, sgh_idx  , 0.0_r8)
      call pbuf_set_field(pbuf2d, sgh30_idx, 0.0_r8)
      call pbuf_set_field(pbuf2d, landm_idx, 0.0_r8)
    end if

    call infld('PBLH', fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, tptr(:,:), found, gridname='physgrid')
    if (.not. found) then
      tptr(:,:) = 0._r8
      if (masterproc) write(iulog, *) 'PBLH initialized to 0.'
    end if
    pblh_idx = pbuf_get_index('pblh')

    call pbuf_set_field(pbuf2d, pblh_idx, tptr)

    call infld('TPERT', fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, tptr(:,:), found, gridname='physgrid')
    if (.not. found) then
      tptr(:,:) = 0._r8
      if (masterproc) write(iulog, *) 'TPERT initialized to 0.'
    end if
    tpert_idx = pbuf_get_index( 'tpert')
    call pbuf_set_field(pbuf2d, tpert_idx, tptr)

    fieldname = 'QPERT'
    qpert_idx = pbuf_get_index( 'qpert',ierr)
    if (qpert_idx > 0) then
      call infld(fieldname, fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, tptr, found, gridname='physgrid')
      if (.not. found) then
        tptr=0_r8
        if (masterproc) write(iulog, *) trim(fieldname), ' initialized to 0.'
      end if

      allocate(tptr3d_2(pcols,pcnst,begchunk:endchunk))
      tptr3d_2 = 0_r8
      tptr3d_2(:,1,:) = tptr(:,:)

      call pbuf_set_field(pbuf2d, qpert_idx, tptr3d_2)
      deallocate(tptr3d_2)
    end if

    fieldname='CUSH'
    m = pbuf_get_index('cush', ierr)
    if (m > 0) then
      call infld(fieldname, fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, tptr, found, gridname='physgrid')
      if (.not.found) then
        if (masterproc) write(iulog, *) trim(fieldname), ' initialized to 1000.'
        tptr = 1000.0_r8
      end if
      do n = 1, dyn_time_lvls
        call pbuf_set_field(pbuf2d, m, tptr, start=[1,n], kount=[pcols,1])
      end do
      deallocate(tptr)
    end if

    !
    ! 3-D fields
    !
    allocate(tptr3d(pcols,pver,begchunk:endchunk))

    fieldname = 'CLOUD'
    m = pbuf_get_index('CLD')
    call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d, found, gridname='physgrid')
    if (found) then
      do n = 1, dyn_time_lvls
        call pbuf_set_field(pbuf2d, m, tptr3d, [1,1,n],[pcols,pver,1])
      end do
    else
      call pbuf_set_field(pbuf2d, m, 0._r8)
      if (masterproc) write(iulog, *) trim(fieldname), ' initialized to 0.'
    end if

    fieldname='QCWAT'
    m = pbuf_get_index(fieldname, ierr)
    if (m > 0) then
      call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d, found, gridname='physgrid')
      if (.not. found) then
        call infld('Q',fh_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d, found, gridname='physgrid')
        if (found) then
          if (masterproc) write(iulog,*) trim(fieldname), ' initialized with Q'
          if(dycore_is('LR')) call polar_average(pver, tptr3d)
        else
          if (masterproc) write(iulog,*) trim(fieldname), ' initialized to huge()'
          tptr3d = huge(1.0_r8)
        end if
      end if
      do n = 1, dyn_time_lvls
        call pbuf_set_field(pbuf2d, m, tptr3d, [1,1,n],[pcols,pver,1])
      end do
    end if

    fieldname = 'ICCWAT'
    m = pbuf_get_index(fieldname, ierr)
    if (m > 0) then
      call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d, found, gridname='physgrid')
      if (found) then
        do n = 1, dyn_time_lvls
          call pbuf_set_field(pbuf2d, m, tptr3d, [1,1,n],[pcols,pver,1])
        end do
      else
        call cnst_get_ind('CLDICE', ixcldice)
        call infld('CLDICE', fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d, found, gridname='physgrid')
        if (found) then
          do n = 1, dyn_time_lvls
            call pbuf_set_field(pbuf2d, m, tptr3d, [1,1,n],[pcols,pver,1])
          end do
        else
          call pbuf_set_field(pbuf2d, m, 0._r8)
        end if
        if (masterproc) then
          if (found) then
            write(iulog, *) trim(fieldname), ' initialized with CLDICE'
          else
            write(iulog, *) trim(fieldname), ' initialized to 0.0'
          end if
        end if
      end if
    end if

    fieldname = 'LCWAT'
    m = pbuf_get_index(fieldname,ierr)
    if (m > 0) then
      call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d, found, gridname='physgrid')
      if (found) then
        do n = 1, dyn_time_lvls
          call pbuf_set_field(pbuf2d, m, tptr3d, [1,1,n],[pcols,pver,1])
        end do
      else
        allocate(tptr3d_2(pcols,pver,begchunk:endchunk))
        call cnst_get_ind('CLDICE', ixcldice)
        call cnst_get_ind('CLDLIQ', ixcldliq)
        call infld('CLDICE', fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d, found, gridname='physgrid')
        call infld('CLDLIQ', fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d_2, found2, gridname='physgrid')
        if(found .and. found2) then
          do lchnk = begchunk, endchunk
            ncol = get_ncols_p(lchnk)
            tptr3d(:ncol,:,lchnk) = tptr3d(:ncol,:,lchnk)+tptr3d_2(:ncol,:,lchnk)
          end do
          if (masterproc) write(iulog, *) trim(fieldname), ' initialized with CLDICE + CLDLIQ'
        else if (found) then ! Data already loaded in tptr3d
          if (masterproc) write(iulog, *) trim(fieldname), ' initialized with CLDICE only'
        else if (found2) then
          tptr3d(:,:,:) = tptr3d_2(:,:,:)
          if (masterproc) write(iulog, *) trim(fieldname), ' initialized with CLDLIQ only'
        end if

        if (found .or. found2) then
          do n = 1, dyn_time_lvls
            call pbuf_set_field(pbuf2d, m, tptr3d, [1,1,n],[pcols,pver,1])
          end do
          if(dycore_is('LR')) call polar_average(pver, tptr3d)
        else
          call pbuf_set_field(pbuf2d, m, 0.0_r8)
          if (masterproc)  write(iulog, *) trim(fieldname), ' initialized to 0.0'
        end if
        deallocate(tptr3d_2)
      end if
    end if

    deallocate(tptr3d)
    allocate(tptr3d(pcols,pver,begchunk:endchunk))

    fieldname = 'TCWAT'
    m = pbuf_get_index(fieldname, ierr)
    if (m > 0) then
      call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d, found, gridname='physgrid')
      if (.not. found) then
        call infld('T', fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d, found, gridname='physgrid')
        if (found) then
          if (dycore_is('LR')) call polar_average(pver, tptr3d)
          if (masterproc) write(iulog, *) trim(fieldname), ' initialized with T'
        else
          if (masterproc) write(iulog, *) trim(fieldname), ' initialized to huge()'
          tptr3d = huge(1._r8)
        end if
      end if
      do n = 1, dyn_time_lvls
        call pbuf_set_field(pbuf2d, m, tptr3d, [1,1,n], [pcols,pver,1])
      end do
    end if

    deallocate(tptr3d)
    allocate(tptr3d(pcols,pverp,begchunk:endchunk))

    fieldname = 'TKE'
    m = pbuf_get_index( 'tke')
    call infld(fieldname, fh_ini, dim1name, 'ilev', dim2name, 1, pcols, 1, pverp, begchunk, endchunk, tptr3d, found, gridname='physgrid')
    if (found) then
      call pbuf_set_field(pbuf2d, m, tptr3d)
    else
      call pbuf_set_field(pbuf2d, m, 0.01_r8)
      if (masterproc) write(iulog, *) trim(fieldname), ' initialized to 0.01'
    end if

    fieldname = 'KVM'
    m = pbuf_get_index('kvm')
    call infld(fieldname, fh_ini, dim1name, 'ilev', dim2name, 1, pcols, 1, pverp, begchunk, endchunk, tptr3d, found, gridname='physgrid')
    if (found) then
      call pbuf_set_field(pbuf2d, m, tptr3d)
    else
      call pbuf_set_field(pbuf2d, m, 0._r8)
      if (masterproc) write(iulog, *) trim(fieldname), ' initialized to 0.'
    end if

    fieldname = 'KVH'
    m = pbuf_get_index('kvh')
    call infld(fieldname, fh_ini, dim1name, 'ilev', dim2name, 1, pcols, 1, pverp, begchunk, endchunk, tptr3d, found, gridname='physgrid')
    if (found) then
      call pbuf_set_field(pbuf2d, m, tptr3d)
    else
      call pbuf_set_field(pbuf2d, m, 0._r8)
      if (masterproc) write(iulog, *) trim(fieldname), ' initialized to 0.'
    end if

    deallocate(tptr3d)
    allocate(tptr3d(pcols,pver,begchunk:endchunk))

    fieldname = 'CONCLD'
    m = pbuf_get_index('CONCLD',ierr)
    if (m > 0) then
      call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, tptr3d, found, gridname='physgrid')
      if(found) then
        do n = 1, dyn_time_lvls
          call pbuf_set_field(pbuf2d, m, tptr3d, [1,1,n],[pcols,pver,1])
        end do
      else
        call pbuf_set_field(pbuf2d, m, 0._r8)
        if (masterproc) write(iulog, *) trim(fieldname), ' initialized to 0.'
      end if

      deallocate (tptr3d)
    end if

    call initialize_short_lived_species(fh_ini, pbuf2d)

    !---------------------------------------------------------------------------------
    !  If needed, get ion and electron temperature fields from initial condition file
    !---------------------------------------------------------------------------------

    call waccmx_phys_ion_elec_temp_inidat(fh_ini, pbuf2d)

  end subroutine phys_inidat

  subroutine phys_init(phys_state, phys_tend, pbuf2d, cam_in, cam_out)

    use physics_buffer      , only: physics_buffer_desc, pbuf_initialize, pbuf_get_index
    use physconst           , only: rair, cpair, gravit, stebol, tmelt, &
                                    latvap, latice, rh2o, rhoh2o, pstd, zvir, &
                                    karman, rhodair, physconst_init
    use ref_pres            , only: pref_edge, pref_mid
    use carma_intr          , only: carma_init
    use cam_control_mod     , only: initial_run
    use check_energy        , only: check_energy_init
    use chemistry           , only: chem_init
    use prescribed_ozone    , only: prescribed_ozone_init
    use prescribed_ghg      , only: prescribed_ghg_init
    use prescribed_aero     , only: prescribed_aero_init
    use aerodep_flx         , only: aerodep_flx_init
    use aircraft_emit       , only: aircraft_emit_init
    use prescribed_volcaero , only: prescribed_volcaero_init
    use prescribed_strataero, only: prescribed_strataero_init
    use cloud_fraction      , only: cldfrc_init
    use cldfrc2m            , only: cldfrc2m_init
    use co2_cycle           , only: co2_init, co2_transport
    use convect_deep        , only: convect_deep_init
    use convect_shallow     , only: convect_shallow_init
    use cam_diagnostics     , only: diag_init
    use gw_drag             , only: gw_init
    use cam3_aero_data      , only: cam3_aero_data_on, cam3_aero_data_init
    use cam3_ozone_data     , only: cam3_ozone_data_on, cam3_ozone_data_init
    use radheat             , only: radheat_init
    use radiation           , only: radiation_init
    use cloud_diagnostics   , only: cloud_diagnostics_init
    use rk_stratiform       , only: rk_stratiform_init
    use wv_saturation       , only: wv_sat_init
    use microp_driver       , only: microp_driver_init
    use microp_aero         , only: microp_aero_init
    use macrop_driver       , only: macrop_driver_init
    use conv_water          , only: conv_water_init
    use aoa_tracers         , only: aoa_tracers_init
    use rayleigh_friction   , only: rayleigh_friction_init
    use pbl_utils           , only: pbl_utils_init
    use vertical_diffusion  , only: vertical_diffusion_init
    use phys_debug_util     , only: phys_debug_init
    use phys_debug          , only: phys_debug_state_init
    use rad_constituents    , only: rad_cnst_init
    use aer_rad_props       , only: aer_rad_props_init
    use subcol              , only: subcol_init
    use qbo                 , only: qbo_init
    use qneg_module         , only: qneg_init
    use lunar_tides         , only: lunar_tides_init
    use iondrag             , only: iondrag_init, do_waccm_ions
    use epp_ionization      , only: epp_ionization_init, epp_ionization_active
    use waccmx_phys_intr    , only: waccmx_phys_ion_elec_temp_init  ! Initialization of ionosphere module (WACCM-X)
    use waccmx_phys_intr    , only: waccmx_phys_mspd_init           ! Initialization of major species diffusion module (WACCM-X)
    use clubb_intr          , only: clubb_ini_cam
    use sslt_rebin          , only: sslt_rebin_init
    use tropopause          , only: tropopause_init
    use solar_data          , only: solar_data_init
    use dadadj_cam          , only: dadadj_init
    use cam_abortutils      , only: endrun

    type(physics_state), pointer       :: phys_state(:)
    type(physics_tend ), pointer       :: phys_tend(:)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)

    type(cam_in_t),  intent(in)        :: cam_in(begchunk:endchunk)
    type(cam_out_t), intent(inout)     :: cam_out(begchunk:endchunk)

    integer lchnk, ierr

    call physics_type_alloc(phys_state, phys_tend, begchunk, endchunk, pcols)

    do lchnk = begchunk, endchunk
      call physics_state_set_grid(lchnk, phys_state(lchnk))
    end do

    ! Initialize any variables in physconst which are not temporally and/or
    ! spatially constant.
    call physconst_init()

    call phys_debug_init()

    call pbuf_initialize(pbuf2d)

    call subcol_init(pbuf2d)

    ! diag_init makes addfld calls for dynamics fields that are output from
    ! the physics decomposition.
    call diag_init(pbuf2d)

    call check_energy_init()

    call aoa_tracers_init()

    teout_idx = pbuf_get_index( 'TEOUT')

    if (adiabatic .or. ideal_phys) then
      if (adiabatic) then
        call endrun('phys_init: adiabatic configuration error')
      else
        call endrun('phys_init: ideal_phys configuration error')
      end if
    end if

    ! if (initial_run) then
    !   call phys_inidat(cam_out, pbuf2d)
    ! end if

    ! wv_saturation is relatively independent of everything else and low level,
    ! so initialize it early. Must at least do this before radiation.
    call wv_sat_init()

    if (cam3_aero_data_on) call cam3_aero_data_init(phys_state)

    call rad_cnst_init()
    call aer_rad_props_init()

    call carma_init()

    call solar_data_init()

    call chem_init(phys_state, pbuf2d)

    call prescribed_ozone_init()
    call prescribed_ghg_init()
    call prescribed_aero_init()
    call aerodep_flx_init()
    call aircraft_emit_init()
    call prescribed_volcaero_init()
    call prescribed_strataero_init()

    if (co2_transport()) call co2_init()

    if (cam3_ozone_data_on) call cam3_ozone_data_init(phys_state)

    call gw_init()

    call rayleigh_friction_init()

    call pbl_utils_init(gravit, karman, cpair, rair, zvir)
    call vertical_diffusion_init(pbuf2d)

    if (waccmx_is('ionosphere') .or. waccmx_is('neutral')) then
      call waccmx_phys_mspd_init ()
      if (waccmx_is('ionosphere')) then
        call waccmx_phys_ion_elec_temp_init(pbuf2d)
      end if
    end if

    call radiation_init(pbuf2d)

    call cloud_diagnostics_init()

    call radheat_init(pref_mid)

    call convect_shallow_init(pref_edge, pbuf2d)

    call cldfrc_init()
    call cldfrc2m_init()

    call convect_deep_init(pref_edge)

    if (microp_scheme == 'RK') then
      call rk_stratiform_init()
    else if (microp_scheme == 'MG') then
      if (.not. do_clubb_sgs) call macrop_driver_init(pbuf2d)
      call microp_aero_init(pbuf2d)
      call microp_driver_init(pbuf2d)
      call conv_water_init()
    else if (microp_scheme == 'SPCAM_m2005') then
      call conv_water_init()
    end if

    if (do_clubb_sgs) call clubb_ini_cam(pbuf2d)

    call qbo_init()

    call lunar_tides_init()

    call iondrag_init(pref_mid)
    if (epp_ionization_active) then
      call epp_ionization_init()
    endif

    call sslt_rebin_init()
    call tropopause_init()
    call dadadj_init()

    prec_dp_idx  = pbuf_get_index('PREC_DP')
    snow_dp_idx  = pbuf_get_index('SNOW_DP')
    prec_sh_idx  = pbuf_get_index('PREC_SH')
    snow_sh_idx  = pbuf_get_index('SNOW_SH')

    dlfzm_idx = pbuf_get_index('DLFZM', ierr)

    call phys_getopts(prog_modal_aero_out=prog_modal_aero)

    if (clim_modal_aero) then
      ! If climate calculations are affected by prescribed modal aerosols, the
      ! the initialization routine for the dry mode radius calculation is called
      ! here.  For prognostic MAM the initialization is called from
      ! modal_aero_initialize
      if (.not. prog_modal_aero) then
         call modal_aero_calcsize_init(pbuf2d)
      end if
      call modal_aero_wateruptake_init(pbuf2d)
    end if

    call qneg_init()

  end subroutine phys_init

  subroutine phys_run1(ztodt, phys_state, phys_tend, pbuf2d,  cam_in, cam_out)

    !-----------------------------------------------------------------------
    !
    ! Purpose:
    ! First part of atmospheric physics package before updating of surface models
    !
    !-----------------------------------------------------------------------
    use time_manager   , only: get_nstep
    use cam_diagnostics, only: diag_allocate, diag_physvar_ic
    use check_energy   , only: check_energy_gmean
    use phys_control   , only: phys_getopts
    use spmd_utils     , only: mpicom
    use physics_buffer , only: physics_buffer_desc, pbuf_get_chunk, pbuf_allocate
#if (defined BFB_CAM_SCAM_IOP)
    use cam_history    , only: outfld
#endif
    use cam_abortutils , only: endrun

    real(r8), intent(in) :: ztodt            ! physics time step unless nstep=0
    type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
    type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend
    type(physics_buffer_desc), pointer, dimension(:,:) :: pbuf2d
    type(cam_in_t),                     dimension(begchunk:endchunk) :: cam_in
    type(cam_out_t),                    dimension(begchunk:endchunk) :: cam_out

    integer c, ncol, nstep
    logical use_spcam
    type(physics_buffer_desc), pointer :: phys_buffer_chunk(:)

    nstep = get_nstep()

    ! The following initialization depends on the import state (cam_in)
    ! being initialized.  This isn't true when cam_init is called, so need
    ! to postpone this initialization to here.
    if (nstep == 0 .and. phys_do_flux_avg()) call flux_avg_init(cam_in,  pbuf2d)

    ! Compute total energy of input state and previous output state
    ! call check_energy_gmean(phys_state, pbuf2d, ztodt, nstep)

    call pbuf_allocate(pbuf2d, 'physpkg')
    call diag_allocate()

    !-----------------------------------------------------------------------
    ! Advance time information
    !-----------------------------------------------------------------------

    call phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d)

#ifdef TRACER_CHECK
    call gmean_mass('before tphysbc DRY', phys_state)
#endif

    !-----------------------------------------------------------------------
    ! Tendency physics before flux coupler invocation
    !-----------------------------------------------------------------------

    call phys_getopts(use_spcam_out=use_spcam)

!$OMP PARALLEL DO PRIVATE (c, phys_buffer_chunk)
    do c = begchunk, endchunk
      !
      ! Output physics terms to IC file
      !
      phys_buffer_chunk => pbuf_get_chunk(pbuf2d, c)

      call diag_physvar_ic(c, phys_buffer_chunk, cam_out(c), cam_in(c))

      call tphysbc(ztodt, phys_state(c), phys_tend(c), phys_buffer_chunk, cam_out(c), cam_in(c))
    end do

#ifdef TRACER_CHECK
    call gmean_mass('between DRY', phys_state)
#endif

  end subroutine phys_run1

  subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d,  cam_out, cam_in)
    !-----------------------------------------------------------------------
    !
    ! Purpose:
    ! Second part of atmospheric physics package after updating of surface models
    !
    !-----------------------------------------------------------------------
    use physics_buffer,  only: physics_buffer_desc, pbuf_get_chunk, pbuf_deallocate, pbuf_update_tim_idx
    use mo_lightning,    only: lightning_no_prod
    use cam_diagnostics, only: diag_deallocate, diag_surf
    use physconst,       only: stebol, latvap
    use carma_intr,      only: carma_accumulate_stats
    use spmd_utils,      only: mpicom

    real(r8), intent(in) :: ztodt                       ! physics time step unless nstep=0
    type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
    type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend
    type(physics_buffer_desc),pointer,  dimension(:,:)               :: pbuf2d
    type(cam_out_t),     intent(inout), dimension(begchunk:endchunk) :: cam_out
    type(cam_in_t),      intent(inout), dimension(begchunk:endchunk) :: cam_in

    integer c, ncol
    type(physics_buffer_desc), pointer, dimension(:) :: phys_buffer_chunk

    !-----------------------------------------------------------------------
    ! Tendency physics after coupler
    ! Not necessary at terminal timestep.
    !-----------------------------------------------------------------------
    ! Set lightning production of NO
    call lightning_no_prod(phys_state, pbuf2d, cam_in)

!$OMP PARALLEL DO PRIVATE (c, NCOL, phys_buffer_chunk)
    do c = begchunk, endchunk
      ncol = get_ncols_p(c)
      phys_buffer_chunk => pbuf_get_chunk(pbuf2d, c)
      !
      ! surface diagnostics for history files
      !
      call diag_surf(cam_in(c), cam_out(c), phys_state(c), phys_buffer_chunk)

      call tphysac(ztodt, cam_in(c), cam_out(c), phys_state(c), phys_tend(c), phys_buffer_chunk)
    end do

#ifdef TRACER_CHECK
    call gmean_mass('after tphysac FV:WET)', phys_state)
#endif

    call carma_accumulate_stats()

    call pbuf_deallocate(pbuf2d, 'physpkg')

    call pbuf_update_tim_idx()
    call diag_deallocate()

  end subroutine phys_run2

  subroutine phys_final(phys_state, phys_tend, pbuf2d)

    use physics_buffer
    use chemistry
    use carma_intr
    use wv_saturation
    use vertical_diffusion
    use aoa_tracers
    use rayleigh_friction
    !-----------------------------------------------------------------------
    !
    ! Purpose:
    ! Finalization of physics package
    !
    !-----------------------------------------------------------------------
    ! Input/output arguments
    type(physics_state), pointer :: phys_state(:)
    type(physics_tend ), pointer :: phys_tend(:)
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)

    if (associated(pbuf2d)) then
      call pbuf_deallocate(pbuf2d, 'global')
      deallocate(pbuf2d)
    end if
    deallocate(phys_state)
    deallocate(phys_tend)
    call chem_final()
    call carma_final()
    call wv_sat_final()
    call vertical_diffusion_final()
    call aoa_tracers_final()
    call rayleigh_friction_final()

  end subroutine phys_final

  subroutine tphysac(ztodt, cam_in, cam_out, state, tend, pbuf)
    !-----------------------------------------------------------------------
    !
    ! Tendency physics after coupling to land, sea, and ice models.
    !
    ! Computes the following:
    !
    !   o Aerosol Emission at Surface
    !   o Source-Sink for Advected Tracers
    !   o Symmetric Turbulence Scheme - Vertical Diffusion
    !   o Rayleigh Friction
    !   o Dry Deposition of Aerosol
    !   o Enforce Charge Neutrality ( Only for WACCM )
    !   o Gravity Wave Drag
    !   o QBO Relaxation ( Only for WACCM )
    !   o Ion Drag ( Only for WACCM )
    !   o Scale Dry Mass Energy
    !-----------------------------------------------------------------------
    use physics_buffer,     only: physics_buffer_desc, pbuf_set_field, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx
    use shr_kind_mod,       only: r8 => shr_kind_r8
    use chemistry,          only: chem_is_active, chem_timestep_tend, chem_emissions
    use cam_diagnostics,    only: diag_phys_tend_writeout
    use gw_drag,            only: gw_tend
    use vertical_diffusion, only: vertical_diffusion_tend
    use rayleigh_friction,  only: rayleigh_friction_tend
    use constituents,       only: cnst_get_ind
    use physics_types,      only: physics_state, physics_tend, physics_ptend, physics_update,    &
                                  physics_dme_adjust, set_dry_to_wet, physics_state_check
    use waccmx_phys_intr,   only: waccmx_phys_mspd_tend  ! WACCM-X major diffusion
    use waccmx_phys_intr,   only: waccmx_phys_ion_elec_temp_tend ! WACCM-X
    use aoa_tracers,        only: aoa_tracers_timestep_tend
    use physconst,          only: rhoh2o, latvap,latice
    use aero_model,         only: aero_model_drydep
    use carma_intr,         only: carma_emission_tend, carma_timestep_tend
    use carma_flags_mod,    only: carma_do_aerosol, carma_do_emission
    use check_energy,       only: check_energy_chng, calc_te_and_aam_budgets
    use check_energy,       only: check_tracers_data, check_tracers_init, check_tracers_chng
    use time_manager,       only: get_nstep
    use cam_abortutils,     only: endrun
    use dycore,             only: dycore_is
    use cam_control_mod,    only: aqua_planet
    use mo_gas_phase_chemdr,only: map2chm
    use clybry_fam,         only: clybry_fam_set
    use charge_neutrality,  only: charge_balance
    use qbo,                only: qbo_relax
    use iondrag,            only: iondrag_calc, do_waccm_ions
    use flux_avg,           only: flux_avg_run
    use unicon_cam,         only: unicon_cam_org_diags
    use cam_history,        only: hist_fld_active
    use qneg_module,        only: qneg4
    use co2_cycle,          only: co2_cycle_set_ptend
    use lunar_tides,        only: lunar_tides_tend

    real(r8), intent(in) :: ztodt
    type(cam_in_t),      intent(inout) :: cam_in
    type(cam_out_t),     intent(inout) :: cam_out
    type(physics_state), intent(inout) :: state
    type(physics_tend ), intent(inout) :: tend
    type(physics_buffer_desc), pointer :: pbuf(:)

    type(check_tracers_data) tracerint ! tracer mass integrals and cummulative boundary fluxes

    type(physics_ptend) ptend          ! indivdual parameterization tendencies
    real(r8) zero(pcols)               ! array of zeros
    integer i, k, m, lchnk, ncol, nstep
    integer yr, mon, day, tod
    integer ixcldice, ixcldliq
    logical labort
    real(r8) tvm       (pcols,pver)       ! virtual temperature
    real(r8) prect     (pcols)            ! total precipitation
    real(r8) surfric   (pcols)            ! surface friction velocity
    real(r8) obklen    (pcols)            ! Obukhov length
    real(r8) fh2o      (pcols)            ! h2o flux to balance source from methane chemistry
    real(r8) flx_heat  (pcols)            ! Heat flux for check_energy_chng.
    real(r8) tmp_q     (pcols,pver)       ! tmp space
    real(r8) tmp_cldliq(pcols,pver)       ! tmp space
    real(r8) tmp_cldice(pcols,pver)       ! tmp space
    real(r8) tmp_trac  (pcols,pver,pcnst) ! tmp space
    real(r8) tmp_pdel  (pcols,pver)       ! tmp space
    real(r8) tmp_ps    (pcols)            ! tmp space

    ! physics buffer fields for total energy and mass adjustment
    integer itim_old, ifld

    real(r8), pointer, dimension(:,:) :: cld
    real(r8), pointer, dimension(:,:) :: qini
    real(r8), pointer, dimension(:,:) :: cldliqini
    real(r8), pointer, dimension(:,:) :: cldiceini
    real(r8), pointer, dimension(:,:) :: dtcore
    real(r8), pointer, dimension(:,:) :: ast     ! relative humidity cloud fraction

    lchnk = state%lchnk
    ncol  = state%ncol
    nstep = get_nstep()

    ! Adjust the surface fluxes to reduce instabilities in near sfc layer
    if (phys_do_flux_avg()) then
      call flux_avg_run(state, cam_in,  pbuf, nstep, ztodt)
    end if

    ! Validate the physics state.
    if (state_debug_checks) call physics_state_check(state, name='before tphysac')

    ! Associate pointers with physics buffer fields
    itim_old = pbuf_old_tim_idx()

    ifld = pbuf_get_index('DTCORE')
    call pbuf_get_field(pbuf, ifld, dtcore, start=[1,1,itim_old], kount=[pcols,pver,1] )
    call pbuf_get_field(pbuf, qini_idx, qini)
    call pbuf_get_field(pbuf, cldliqini_idx, cldliqini)
    call pbuf_get_field(pbuf, cldiceini_idx, cldiceini)

    ifld = pbuf_get_index('CLD')
    call pbuf_get_field(pbuf, ifld, cld, start=[1,1,itim_old], kount=[pcols,pver,1])

    ifld = pbuf_get_index('AST')
    call pbuf_get_field(pbuf, ifld, ast, start=[1,1,itim_old], kount=[pcols,pver,1])

    !
    ! Accumulate fluxes into net flux array for spectral dycores
    ! jrm Include latent heat of fusion for snow
    !
    do i = 1, ncol
      tend%flx_net(i) = tend%flx_net(i) + cam_in%shf(i) + (cam_out%precc(i) &
                      + cam_out%precl(i)) * latvap * rhoh2o &
                      + (cam_out%precsc(i) + cam_out%precsl(i)) * latice * rhoh2o
    end do

    ! Emissions of aerosols and gas-phase chemistry constituents at surface
    call chem_emissions(state, cam_in)

    if (carma_do_emission) then
      ! carma emissions
      call carma_emission_tend (state, ptend, cam_in, ztodt)
      call physics_update(state, ptend, ztodt, tend)
    end if

    ! get nstep and zero array for energy checker
    zero = 0
    nstep = get_nstep()
    call tracerint%init()
    call check_tracers_init(state, tracerint)

    ! Check if latent heat flux exceeds the total moisture content of the
    ! lowest model layer, thereby creating negative moisture.

    call qneg4('TPHYSAC', lchnk, ncol, ztodt ,   &
         state%q(1,pver,1), state%rpdel(1,pver), &
         cam_in%shf, cam_in%lhf, cam_in%cflx)

    ! --------------------------------------------------------------------------
    ! Source/sink terms for advected tracers.
    ! --------------------------------------------------------------------------

    call aoa_tracers_timestep_tend(state, ptend, cam_in%cflx, cam_in%landfrac, ztodt)
    call physics_update(state, ptend, ztodt, tend)
    call check_tracers_chng(state, tracerint, 'aoa_tracers_timestep_tend', nstep, ztodt, cam_in%cflx)

    call co2_cycle_set_ptend(state, pbuf, ptend)
    call physics_update(state, ptend, ztodt, tend)

    ! --------------------------------------------------------------------------
    ! Chemistry and MAM calculation
    ! MAM core aerosol conversion process is performed in the below 'chem_timestep_tend'.
    ! In addition, surface flux of aerosol species other than 'dust' and 'sea salt', and
    ! elevated emission of aerosol species are treated in 'chem_timestep_tend' before
    ! Gas chemistry and MAM core aerosol conversion.
    ! Note that surface flux is not added into the atmosphere, but elevated emission is
    ! added into the atmosphere as tendency.
    ! --------------------------------------------------------------------------
    if (chem_is_active()) then
      call chem_timestep_tend(state, ptend, cam_in, cam_out, ztodt, pbuf, fh2o=fh2o)
      call physics_update(state, ptend, ztodt, tend)

      call check_energy_chng(state, tend, 'chem', nstep, ztodt, fh2o, zero, zero, zero)
      call check_tracers_chng(state, tracerint, 'chem_timestep_tend', nstep, ztodt, cam_in%cflx)
    end if

    ! --------------------------------------------------------------------------
    ! Vertical diffusion calculation
    ! Call vertical diffusion code (PBL, free atmosphere and molecular)
    ! --------------------------------------------------------------------------
    call vertical_diffusion_tend(ztodt ,state , cam_in, surfric, obklen, ptend, ast, pbuf)

    ! --------------------------------------------------------------------------
    ! Call major diffusion for extended model
    ! -------------------------------------------------------------------------- 
    if (waccmx_is('ionosphere') .or. waccmx_is('neutral')) then
      call waccmx_phys_mspd_tend(ztodt, state, ptend)
    end if

    call physics_update(state, ptend, ztodt, tend)

    ! --------------------------------------------------------------------------
    ! Rayleigh friction calculation
    ! --------------------------------------------------------------------------
    call rayleigh_friction_tend(ztodt, state, ptend)
    call physics_update(state, ptend, ztodt, tend)

    if (do_clubb_sgs) then
      call check_energy_chng(state, tend, 'vdiff', nstep, ztodt, zero, zero, zero, zero)
    else
      call check_energy_chng(state, tend, 'vdiff', nstep, ztodt, cam_in%cflx(:,1), zero, zero, cam_in%shf)
    end if

    call check_tracers_chng(state, tracerint, 'vdiff', nstep, ztodt, cam_in%cflx)

    ! Aerosol dry deposition processes

    call aero_model_drydep(state, pbuf, obklen, surfric, cam_in, ztodt, cam_out, ptend)
    call physics_update(state, ptend, ztodt, tend)

    ! CARMA microphysics
    !
    ! NOTE: This does both the timestep_tend for CARMA aerosols as well as doing the dry
    ! deposition for CARMA aerosols. It needs to follow vertical_diffusion_tend, so that
    ! obklen and surfric have been calculated. It needs to follow aero_model_drydep, so
    ! that cam_out%xxxdryxxx fields have already been set for CAM aerosols and cam_out
    ! can be added to for CARMA aerosols.
    if (carma_do_aerosol) then
      call carma_timestep_tend(state, cam_in, cam_out, ptend, ztodt, pbuf, obklen=obklen, ustar=surfric)
      call physics_update(state, ptend, ztodt, tend)
      call check_energy_chng(state, tend, 'carma_tend', nstep, ztodt, zero, zero, zero, zero)
    end if

    ! --------------------------------------------------------------------------
    ! Enforce charge neutrality
    ! --------------------------------------------------------------------------
    call charge_balance(state, pbuf)

    ! --------------------------------------------------------------------------
    ! Gravity wave drag
    ! --------------------------------------------------------------------------

    call gw_tend(state, pbuf, ztodt, ptend, cam_in, flx_heat)
    call physics_update(state, ptend, ztodt, tend)

    ! Check energy integrals
    call check_energy_chng(state, tend, 'gwdrag', nstep, ztodt, zero, zero, zero, flx_heat)

    ! --------------------------------------------------------------------------
    ! QBO relaxation
    ! --------------------------------------------------------------------------

    call qbo_relax(state, pbuf, ptend)
    call physics_update(state, ptend, ztodt, tend)

    ! Check energy integrals
    call check_energy_chng(state, tend, 'qborelax', nstep, ztodt, zero, zero, zero, zero)

    ! --------------------------------------------------------------------------
    ! Lunar tides
    ! --------------------------------------------------------------------------

    call lunar_tides_tend( state, ptend )
    call physics_update(state, ptend, ztodt, tend)

    ! Check energy integrals
    call check_energy_chng(state, tend, 'lunar_tides', nstep, ztodt, zero, zero, zero, zero)

    ! --------------------------------------------------------------------------
    ! Ion drag calculation
    ! --------------------------------------------------------------------------

    if (do_waccm_ions) then
      call iondrag_calc(lchnk, ncol, state, ptend, pbuf, ztodt)
    else
      call iondrag_calc(lchnk, ncol, state, ptend)
    end if

    !----------------------------------------------------------------------------
    ! Call ionosphere routines for extended model if mode is set to ionosphere
    !----------------------------------------------------------------------------
    if (waccmx_is('ionosphere')) then
      call waccmx_phys_ion_elec_temp_tend(state, ptend, pbuf, ztodt)
    end if

    call physics_update(state, ptend, ztodt, tend)

    call calc_te_and_aam_budgets(state, 'pAP')

    !---------------------------------------------------------------------------------
    ! Enforce charge neutrality after O+ change from ionos_tend
    !---------------------------------------------------------------------------------
    if (waccmx_is('ionosphere')) then
      call charge_balance(state, pbuf)
    end if

    ! Check energy integrals
    call check_energy_chng(state, tend, 'iondrag', nstep, ztodt, zero, zero, zero, zero)

    ! Energy budget checks

    ! Save total energy for global fixer in next timestep (FV and SE dycores)
    call pbuf_set_field(pbuf, teout_idx, state%te_cur, [1,itim_old],[pcols,1])

    if (shallow_scheme == 'UNICON') then

      ! ------------------------------------------------------------------------
      ! Insert the organization-related heterogeneities computed inside the
      ! UNICON into the tracer arrays here before performing advection.
      ! This is necessary to prevent any modifications of organization-related
      ! heterogeneities by non convection-advection process, such as
      ! dry and wet deposition of aerosols, MAM, etc.
      ! Again, note that only UNICON and advection schemes are allowed to
      ! changes to organization at this stage, although we can include the
      ! effects of other physical processes in future.
      ! ------------------------------------------------------------------------

      call unicon_cam_org_diags(state, pbuf)

    end if
    !
    ! FV: convert dry-type mixing ratios to moist here because physics_dme_adjust
    !     assumes moist. This is done in p_d_coupling for other dynamics. Bundy, Feb 2004.
    if (dycore_is('LR') .or. dycore_is('FV3')) call set_dry_to_wet(state)    ! Physics had dry, dynamics wants moist

    ! Scale dry mass and energy (does nothing if dycore is EUL or SLD)
    call cnst_get_ind('CLDLIQ', ixcldliq)
    call cnst_get_ind('CLDICE', ixcldice)

    tmp_q     (:ncol,:pver) = state%q(:ncol,:pver,1)
    tmp_cldliq(:ncol,:pver) = state%q(:ncol,:pver,ixcldliq)
    tmp_cldice(:ncol,:pver) = state%q(:ncol,:pver,ixcldice)

    ! For not ('FV'|'FV3'), physics_dme_adjust is called for energy diagnostic purposes only.  So, save off tracers
    if (.not. (dycore_is('FV') .or. dycore_is('FV3')) .and. &
        (hist_fld_active('SE_pAM') .or. hist_fld_active('KE_pAM') .or. hist_fld_active('WV_pAM') .or. &
         hist_fld_active('WL_pAM') .or. hist_fld_active('WI_pAM'))) then
      tmp_trac(:ncol,:pver,:pcnst) = state%q   (:ncol,:pver,:pcnst)
      tmp_pdel(:ncol,:pver)        = state%pdel(:ncol,:pver)
      tmp_ps  (:ncol)              = state%ps  (:ncol)
      !
      ! pint, lnpint,rpdel are altered by dme_adjust but not used for tendencies in dynamics of SE
      ! we do not reset them to pre-dme_adjust values
      !
      if (dycore_is('SE')) call set_dry_to_wet(state)

      call physics_dme_adjust(state, tend, qini, ztodt)

      call calc_te_and_aam_budgets(state, 'pAM')
      ! Restore pre-'physics_dme_adjust' tracers
      state%q(:ncol,:pver,:pcnst) = tmp_trac(:ncol,:pver,:pcnst)
      state%pdel(:ncol,:pver)     = tmp_pdel(:ncol,:pver)
      state%ps(:ncol)             = tmp_ps(:ncol)
    end if

    if (dycore_is('LR') .or. dycore_is('FV3')) then
      call physics_dme_adjust(state, tend, qini, ztodt)
      call calc_te_and_aam_budgets(state, 'pAM')
    end if

    ! Store T in buffer for use in computing dynamics T-tendency in next timestep
    do k = 1, pver
      dtcore(:ncol,k) = state%t(:ncol,k)
    end do

    if (aqua_planet) then
      labort = .false.
      do i = 1, ncol
        if (cam_in%ocnfrac(i) /= 1.0_r8) then
          labort = .true.
          if (masterproc) write(iulog,*) 'oceanfrac(',i,')=',cam_in%ocnfrac(i)
        end if
      end do
      if (labort) then
        call endrun('TPHYSAC error: in aquaplanet mode, but grid contains non-ocean point')
      end if
    end if

    call diag_phys_tend_writeout(state, pbuf,  tend, ztodt, tmp_q, tmp_cldliq, tmp_cldice, &
                                 qini, cldliqini, cldiceini)

    call clybry_fam_set(ncol, lchnk, map2chm, state%q, pbuf)

  end subroutine tphysac

  subroutine tphysbc(ztodt, state, tend, pbuf, cam_out, cam_in)

    !-----------------------------------------------------------------------
    !
    ! Purpose:
    !
    !   Evaluate and apply physical processes that are calculated before
    !   coupling to land, sea, and ice models.
    !
    ! Processes currently included are:
    !
    !   o Resetting Negative Tracers to Positive
    !   o Global Mean Total Energy Fixer
    !   o Dry Adjustment
    !   o Asymmetric Turbulence Scheme : Deep Convection & Shallow Convection
    !   o Stratiform Macro-Microphysics
    !   o Wet Scavenging of Aerosol
    !   o Radiation
    !
    ! Method:
    !
    !   Each parameterization should be implemented with this sequence of calls:
    !    1)  Call physics interface
    !    2)  Check energy
    !    3)  Call physics_update
    !   See Interface to Column Physics and Chemistry Packages
    !     http://www.ccsm.ucar.edu/models/atm-cam/docs/phys-interface/index.html
    !
    !-----------------------------------------------------------------------

    use physics_buffer     , only: physics_buffer_desc, pbuf_get_field
    use physics_buffer     , only: pbuf_get_index, pbuf_old_tim_idx
    use physics_buffer     , only: col_type_subcol, dyn_time_lvls
    use shr_kind_mod       , only: r8 => shr_kind_r8
    use dadadj_cam         , only: dadadj_tend
    use rk_stratiform      , only: rk_stratiform_tend
    use microp_driver      , only: microp_driver_tend
    use microp_aero        , only: microp_aero_run
    use macrop_driver      , only: macrop_driver_tend
    use physics_types      , only: physics_state, physics_tend, physics_ptend, &
                                   physics_update, physics_ptend_init, physics_ptend_sum, &
                                   physics_state_check, physics_ptend_scale
    use cam_diagnostics    , only: diag_conv_tend_ini, diag_phys_writeout, diag_conv, diag_export, diag_state_b4_phys_write
    use cam_diagnostics    , only: diag_clip_tend_writeout
    use cam_history        , only: outfld
    use physconst          , only: cpair, latvap
    use constituents       , only: pcnst, qmin, cnst_get_ind
    use convect_deep       , only: convect_deep_tend, convect_deep_tend_2, deep_scheme_does_scav_trans
    use time_manager       , only: is_first_step, get_nstep
    use convect_shallow    , only: convect_shallow_tend
    use check_energy       , only: check_energy_chng, check_energy_fix, check_energy_timestep_init
    use check_energy       , only: check_tracers_data, check_tracers_init, check_tracers_chng
    use check_energy       , only: calc_te_and_aam_budgets
    use dycore             , only: dycore_is
    use aero_model         , only: aero_model_wetdep
    use carma_intr         , only: carma_wetdep_tend, carma_timestep_tend
    use carma_flags_mod    , only: carma_do_detrain, carma_do_cldice, carma_do_cldliq,  carma_do_wetdep
    use radiation          , only: radiation_tend
    use cloud_diagnostics  , only: cloud_diagnostics_calc
    use mo_gas_phase_chemdr, only: map2chm
    use clubb_intr         , only: clubb_tend_cam
    use sslt_rebin         , only: sslt_rebin_adv
    use tropopause         , only: tropopause_output
    use cam_abortutils     , only: endrun
    use subcol             , only: subcol_gen, subcol_ptend_avg
    use subcol_utils       , only: subcol_ptend_copy, is_subcol_on
    use qneg_module        , only: qneg3
    use subcol_SILHS       , only: subcol_SILHS_var_covar_driver
    use subcol_SILHS       , only: subcol_SILHS_fill_holes_conserv
    use subcol_SILHS       , only: subcol_SILHS_hydromet_conc_tend_lim
    use micro_mg_cam       , only: massless_droplet_destroyer

    real(r8), intent(in) :: ztodt ! 2 delta t (model time increment)

    type(physics_state), intent(inout) :: state
    type(physics_tend ), intent(inout) :: tend
    type(physics_buffer_desc), pointer :: pbuf(:)
    type(cam_out_t), intent(inout) :: cam_out
    type(cam_in_t), intent(in) :: cam_in

    type(physics_ptend) ptend            ! indivdual parameterization tendencies
    type(physics_state) state_sc         ! state for sub-columns
    type(physics_ptend) ptend_sc         ! ptend for sub-columns
    type(physics_ptend) ptend_aero       ! ptend for microp_aero
    type(physics_ptend) ptend_aero_sc    ! ptend for microp_aero on sub-columns
    type(physics_tend)  tend_sc          ! tend for sub-columns

    real(r8) net_flx(pcols)
    real(r8) zdu(pcols,pver)             ! detraining mass flux from deep convection
    real(r8) cmfmc(pcols,pverp)          ! Convective mass flux--m sub c
    real(r8) cmfcme(pcols,pver)          ! cmf condensation - evaporation
    real(r8) dlf(pcols,pver)             ! Detraining cld H20 from shallow + deep convections
    real(r8) dlf2(pcols,pver)            ! Detraining cld H20 from shallow convections
    real(r8) pflx(pcols,pverp)           ! Conv rain flux thru out btm of lev
    real(r8) rtdt                        ! 1./ztodt
    integer lchnk, ncol, nstep, i
    integer ixcldice, ixcldliq, ixq      ! constituent indices for cloud liquid and ice water.
    ! for macro/micro co-substepping
    integer macmic_it                    ! iteration variables
    real(r8) cld_macmic_ztodt            ! modified timestep
    ! physics buffer fields to compute tendencies for stratiform package
    integer itim_old, ifld
    real(r8), pointer, dimension(:,:) :: cld       ! cloud fraction

    ! physics buffer fields for total energy and mass adjustment
    real(r8), pointer, dimension(:  ) :: teout
    real(r8), pointer, dimension(:,:) :: qini
    real(r8), pointer, dimension(:,:) :: cldliqini
    real(r8), pointer, dimension(:,:) :: cldiceini
    real(r8), pointer, dimension(:,:) :: dtcore

    real(r8), pointer, dimension(:,:,:) :: fracis  ! fraction of transported species that are insoluble

    real(r8), pointer :: dlfzm(:,:)      ! ZM detrained convective cloud water mixing ratio.

    ! convective precipitation variables
    real(r8), pointer :: prec_dp(:)      ! total precipitation from ZM convection
    real(r8), pointer :: snow_dp(:)      ! snow from ZM convection
    real(r8), pointer :: prec_sh(:)      ! total precipitation from Hack convection
    real(r8), pointer :: snow_sh(:)      ! snow from Hack convection

    ! carma precipitation variables
    real(r8) :: prec_sed_carma(pcols)    ! total precip from cloud sedimentation (CARMA)
    real(r8) :: snow_sed_carma(pcols)    ! snow from cloud ice sedimentation (CARMA)

    ! stratiform precipitation variables
    real(r8), pointer :: prec_str(:)     ! sfc flux of precip from stratiform (m/s)
    real(r8), pointer :: snow_str(:)     ! sfc flux of snow from stratiform   (m/s)
    real(r8), pointer :: prec_str_sc(:)  ! sfc flux of precip from stratiform (m/s) -- for subcolumns
    real(r8), pointer :: snow_str_sc(:)  ! sfc flux of snow from stratiform   (m/s) -- for subcolumns
    real(r8), pointer :: prec_pcw(:)     ! total precip from prognostic cloud scheme
    real(r8), pointer :: snow_pcw(:)     ! snow from prognostic cloud scheme
    real(r8), pointer :: prec_sed(:)     ! total precip from cloud sedimentation
    real(r8), pointer :: snow_sed(:)     ! snow from cloud ice sedimentation

    ! Local copies for substepping
    real(r8) prec_pcw_macmic(pcols)
    real(r8) snow_pcw_macmic(pcols)
    real(r8) prec_sed_macmic(pcols)
    real(r8) snow_sed_macmic(pcols)

    ! energy checking variables
    real(r8) zero    (pcols)             ! array of zeros
    real(r8) zero_sc (pcols*psubcols)    ! array of zeros
    real(r8) rliq    (pcols)             ! vertical integral of liquid not yet in q(ixcldliq)
    real(r8) rice    (pcols)             ! vertical integral of ice not yet in q(ixcldice)
    real(r8) rliq2   (pcols)             ! vertical integral of liquid from shallow scheme
    real(r8) det_s   (pcols)             ! vertical integral of detrained static energy from ice
    real(r8) det_ice (pcols)             ! vertical integral of detrained ice
    real(r8) flx_cnd (pcols)
    real(r8) flx_heat(pcols)
    type(check_tracers_data) tracerint   ! energy integrals and cummulative boundary fluxes
    real(r8) zero_tracers(pcols,pcnst)

    logical lq(pcnst)

    zero         = 0
    zero_tracers = 0
    zero_sc      = 0

    lchnk = state%lchnk
    ncol  = state%ncol
    rtdt  = 1.0_r8 / ztodt
    nstep = get_nstep()

    ! Associate pointers with physics buffer fields
    itim_old = pbuf_old_tim_idx()
    ifld = pbuf_get_index('CLD')
    call pbuf_get_field(pbuf, ifld, cld, [1,1,itim_old], [pcols,pver,1])
    call pbuf_get_field(pbuf, teout_idx, teout, [1,itim_old], [pcols,1])
    call pbuf_get_field(pbuf, qini_idx, qini)
    call pbuf_get_field(pbuf, cldliqini_idx, cldliqini)
    call pbuf_get_field(pbuf, cldiceini_idx, cldiceini)

    ifld = pbuf_get_index('DTCORE')
    call pbuf_get_field(pbuf, ifld, dtcore, start=[1,1,itim_old], kount=[pcols,pver,1])

    ifld = pbuf_get_index('FRACIS')
    call pbuf_get_field(pbuf, ifld, fracis, start=[1,1,1], kount=[pcols, pver, pcnst])
    fracis(:ncol,:,1:pcnst) = 1

    ! Set physics tendencies to 0
    tend%dTdt(:ncol,:pver) = 0
    tend%dudt(:ncol,:pver) = 0
    tend%dvdt(:ncol,:pver) = 0

    ! Verify state coming from the dynamics
    if (state_debug_checks) call physics_state_check(state, name='before tphysbc (dycore?)')

    ! Compute mass integrals of input tracers state
    call check_tracers_init(state, tracerint)

    ! --------------------------------------------------------------------------
    ! Global mean total energy fixer
    ! --------------------------------------------------------------------------

    call calc_te_and_aam_budgets(state, 'pBF')
    if (dycore_is('LR') .or. dycore_is('FV3') .or. dycore_is('SE'))  then
      call check_energy_fix(state, ptend, nstep, flx_heat)
      call physics_update(state, ptend, ztodt, tend)
      call check_energy_chng(state, tend, 'chkengyfix', nstep, ztodt, zero, zero, zero, flx_heat)
      call outfld('EFIX', flx_heat, pcols, lchnk)
    end if
    call calc_te_and_aam_budgets(state, 'pBP')

    ! Save state for convective tendency calculations.
    call diag_conv_tend_ini(state, pbuf)

    call cnst_get_ind('Q', ixq)
    call cnst_get_ind('CLDLIQ', ixcldliq)
    call cnst_get_ind('CLDICE', ixcldice)
    qini     (:ncol,:pver) = state%q(:ncol,:pver,       1)
    cldliqini(:ncol,:pver) = state%q(:ncol,:pver,ixcldliq)
    cldiceini(:ncol,:pver) = state%q(:ncol,:pver,ixcldice)

    call outfld('TEOUT', teout       , pcols, lchnk)
    call outfld('TEINP', state%te_ini, pcols, lchnk)
    call outfld('TEFIX', state%te_cur, pcols, lchnk)

    ! T tendency due to dynamics
    if (nstep > dyn_time_lvls - 1) then
      dtcore(:ncol,:pver) = (state%t(:ncol,:pver) - dtcore(:ncol,:pver)) / ztodt
      call outfld('DTCORE', dtcore, pcols, lchnk)
    end if

    ! --------------------------------------------------------------------------
    ! Dry adjustment
    ! --------------------------------------------------------------------------

    call dadadj_tend(ztodt, state, ptend)
    call physics_update(state, ptend, ztodt, tend)

    ! --------------------------------------------------------------------------
    ! Moist convection
    ! --------------------------------------------------------------------------

    call convect_deep_tend(     &
      mcon    =cmfmc          , &
      cme     =cmfcme         , &
      pflx    =pflx           , &
      zdu     =zdu            , &
      rliq    =rliq           , &
      rice    =rice           , &
      ztodt   =ztodt          , &
      state   =state          , &
      ptend   =ptend          , &
      landfrac=cam_in%landfrac, &
      pbuf    =pbuf)
    call physics_update(state, ptend, ztodt, tend)

    call pbuf_get_field(pbuf, prec_dp_idx , prec_dp )
    call pbuf_get_field(pbuf, snow_dp_idx , snow_dp )
    call pbuf_get_field(pbuf, prec_sh_idx , prec_sh )
    call pbuf_get_field(pbuf, snow_sh_idx , snow_sh )
    call pbuf_get_field(pbuf, prec_str_idx, prec_str)
    call pbuf_get_field(pbuf, snow_str_idx, snow_str)
    call pbuf_get_field(pbuf, prec_sed_idx, prec_sed)
    call pbuf_get_field(pbuf, snow_sed_idx, snow_sed)
    call pbuf_get_field(pbuf, prec_pcw_idx, prec_pcw)
    call pbuf_get_field(pbuf, snow_pcw_idx, snow_pcw)

    if (use_subcol_microp) then
      call pbuf_get_field(pbuf, prec_str_idx, prec_str_sc, col_type=col_type_subcol)
      call pbuf_get_field(pbuf, snow_str_idx, snow_str_sc, col_type=col_type_subcol)
    end if

    ! Check energy integrals, including 'reserved liquid'
    flx_cnd(:ncol) = prec_dp(:ncol) + rliq(:ncol)
    snow_dp(:ncol) = snow_dp(:ncol) + rice(:ncol)
    call check_energy_chng(state, tend, 'convect_deep', nstep, ztodt, zero, flx_cnd, snow_dp, zero)
    snow_dp(:ncol) = snow_dp(:ncol) - rice(:ncol)

    !
    ! Call Hack (1994) convection scheme to deal with shallow/mid-level convection
    !
    if (dlfzm_idx > 0) then
      call pbuf_get_field(pbuf, dlfzm_idx, dlfzm)
      dlf(:ncol,:) = dlfzm(:ncol,:)
    else
      dlf = 0
    end if

    call convect_shallow_tend( &
      ztodt     =ztodt       , &
      cmfmc     =cmfmc       , &
      qc        =dlf         , &
      qc2       =dlf2        , &
      rliq      =rliq        , &
      rliq2     =rliq2       , &
      state     =state       , &
      ptend_all =ptend       , &
      pbuf      =pbuf        , &
      cam_in    =cam_in      )
    call physics_update(state, ptend, ztodt, tend)

    flx_cnd(:ncol) = prec_sh(:ncol) + rliq2(:ncol)
    call check_energy_chng(state, tend, 'convect_shallow', nstep, ztodt, zero, flx_cnd, snow_sh, zero)

    call check_tracers_chng(state, tracerint, 'convect_shallow', nstep, ztodt, zero_tracers)

    ! Rebin the 4-bin version of sea salt into bins for coarse and accumulation
    ! modes that correspond to the available optics data.  This is only necessary
    ! for CAM-RT.  But it's done here so that the microphysics code which is called
    ! from the stratiform interface has access to the same aerosols as the radiation
    ! code.
    call sslt_rebin_adv(pbuf, state)

    !===================================================
    ! Calculate tendencies from CARMA bin microphysics.
    !===================================================
    !
    ! If CARMA is doing detrainment, then on output, rliq no longer represents water reserved
    ! for detrainment, but instead represents potential snow fall. The mass and number of the
    ! snow are stored in the physics buffer and will be incorporated by the MG microphysics.
    !
    ! Currently CARMA cloud microphysics is only supported with the MG microphysics.

    if (carma_do_cldice .or. carma_do_cldliq) then
      call carma_timestep_tend(state, cam_in, cam_out, ptend, ztodt, pbuf, dlf=dlf, rliq=rliq, &
        prec_str=prec_str, snow_str=snow_str, prec_sed=prec_sed_carma, snow_sed=snow_sed_carma)
      call physics_update(state, ptend, ztodt, tend)

      ! Before the detrainment, the reserved condensate is all liquid, but if CARMA is doing
      ! detrainment, then the reserved condensate is snow.
      if (carma_do_detrain) then
        call check_energy_chng(state, tend, 'carma_tend', nstep, ztodt, zero, prec_str+rliq, snow_str+rliq, zero)
      else
        call check_energy_chng(state, tend, 'carma_tend', nstep, ztodt, zero, prec_str, snow_str, zero)
      end if
    end if

    if (microp_scheme == 'RK') then

      !===================================================
      ! Calculate stratiform tendency (sedimentation, detrain, cloud fraction and microphysics )
      !===================================================

      call rk_stratiform_tend(    &
        state=state             , &
        ptend_all=ptend         , &
        pbuf=pbuf               , &
        dtime=ztodt             , &
        icefrac=cam_in%icefrac  , &
        landfrac=cam_in%landfrac, &
        ocnfrac=cam_in%ocnfrac  , &
        snowh=cam_in%snowhland  , &
        dlf=dlf                 , &
        dlf2=dlf2               , &
        rliq=rliq               , &
        cmfmc=cmfmc             , &
        ts=cam_in%ts            , &
        sst=cam_in%sst          , &
        zdu=zdu                 )

       call physics_update(state, ptend, ztodt, tend)
       call check_energy_chng(state, tend, 'cldwat_tend', nstep, ztodt, zero, prec_str, snow_str, zero)

    else if (microp_scheme == 'MG') then
      ! Start co-substepping of macrophysics and microphysics
      cld_macmic_ztodt = ztodt / cld_macmic_num_steps

      ! Clear precip fields that should accumulate.
      prec_sed_macmic = 0
      snow_sed_macmic = 0
      prec_pcw_macmic = 0
      snow_pcw_macmic = 0

      do macmic_it = 1, cld_macmic_num_steps

        !===================================================
        ! Calculate macrophysical tendency (sedimentation, detrain, cloud fraction)
        !===================================================

        ! Don't call Park macrophysics if CLUBB is called
        if (macrop_scheme /= 'CLUBB_SGS') then
          call macrop_driver_tend(     &
            state   =state           , &
            ptend   =ptend           , &
            dtime   =cld_macmic_ztodt, &
            landfrac=cam_in%landfrac , &
            ocnfrac =cam_in%ocnfrac  , &
            snowh   =cam_in%snowhland, &
            dlf     =dlf             , &
            dlf2    =dlf2            , &
            cmfmc   =cmfmc           , &
            ts      =cam_in%ts       , &
            sst     =cam_in%sst      , &
            zdu     =zdu             , &
            pbuf    =pbuf            , &
            det_s   =det_s           , &
            det_ice =det_ice         )

          ! Since we 'added' the reserved liquid back in this routine, we need
          ! to account for it in the energy checker
          flx_cnd (:ncol) = -1.0_r8 * rliq(:ncol)
          flx_heat(:ncol) = det_s(:ncol)

          ! Unfortunately, physics_update does not know what time period
          ! 'tend' is supposed to cover, and therefore can't update it
          ! with substeps correctly. For now, work around this by scaling
          ! ptend down by the number of substeps, then applying it for
          ! the full time (ztodt).
          call physics_ptend_scale(ptend, 1.0_r8 / cld_macmic_num_steps, ncol)
          call physics_update(state, ptend, ztodt, tend)

          call check_energy_chng(state, tend, 'macrop_tend', nstep, ztodt, &
            zero, flx_cnd(:ncol)/cld_macmic_num_steps, &
            det_ice(:ncol)/cld_macmic_num_steps, &
            flx_heat(:ncol)/cld_macmic_num_steps)

        else ! Calculate CLUBB macrophysics

          ! =====================================================
          !    CLUBB call (PBL, shallow convection, macrophysics)
          ! =====================================================

          call clubb_tend_cam(                         &
            state               =state               , &
            ptend_all           =ptend               , &
            pbuf                =pbuf                , &
            hdtime              =cld_macmic_ztodt    , &
            cmfmc               =cmfmc               , &
            cam_in              =cam_in              , &
            macmic_it           =macmic_it           , &
            cld_macmic_num_steps=cld_macmic_num_steps, &
            dlf                 =dlf                 , &
            det_s               =det_s               , &
            det_ice             =det_ice             )

          ! Since we 'added' the reserved liquid back in this routine, we need
          ! to account for it in the energy checker
          flx_cnd (:ncol) = -1.0_r8 * rliq(:ncol)
          flx_heat(:ncol) = cam_in%shf(:ncol) + det_s(:ncol)

          ! Unfortunately, physics_update does not know what time period
          ! 'tend' is supposed to cover, and therefore can't update it
          ! with substeps correctly. For now, work around this by scaling
          ! ptend down by the number of substeps, then applying it for
          ! the full time (ztodt).
          call physics_ptend_scale(ptend, 1.0_r8 / cld_macmic_num_steps, ncol)

          ! Update physics tendencies and copy state to state_eq, because that is
          ! input for microphysics
          call physics_update(state, ptend, ztodt, tend)

          ! Use actual qflux (not lhf/latvap) for consistency with surface fluxes and revised code
          call check_energy_chng(state, tend, 'clubb_tend', nstep, ztodt, &
            cam_in%cflx(:ncol,1) / cld_macmic_num_steps, &
            flx_cnd    (:ncol)   / cld_macmic_num_steps, &
            det_ice    (:ncol)   / cld_macmic_num_steps, &
            flx_heat   (:ncol)   / cld_macmic_num_steps)

        end if

        !===================================================
        ! Calculate cloud microphysics
        !===================================================

        if (is_subcol_on()) then
          ! Allocate sub-column structures.
          call physics_state_alloc(state_sc, lchnk, psubcols*pcols)
          call physics_tend_alloc(tend_sc, psubcols*pcols)

          ! Generate sub-columns using the requested scheme
          call subcol_gen(state, tend, state_sc, tend_sc, pbuf)

          !Initialize check energy for subcolumns
          call check_energy_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol)
        end if

        call microp_aero_run(state, ptend_aero, cld_macmic_ztodt, pbuf)

        if (use_subcol_microp) then

          call microp_driver_tend(state_sc, ptend_sc, cld_macmic_ztodt, pbuf)
          ! Parameterize subcolumn effects on covariances, if enabled
          if (trim(subcol_scheme) == 'SILHS') &
            call subcol_SILHS_var_covar_driver( cld_macmic_ztodt, state_sc, ptend_sc, pbuf )

          ! Average the sub-column ptend for use in gridded update - will not contain ptend_aero
          call subcol_ptend_avg(ptend_sc, state_sc%ngrdcol, lchnk, ptend)

          ! Call the conservative hole filler.
          ! Hole filling is only necessary when using subcolumns.
          ! Note:  this needs to be called after subcol_ptend_avg but before
          !        physics_ptend_scale.
          if (trim(subcol_scheme) == 'SILHS') &
            call subcol_SILHS_fill_holes_conserv(state, cld_macmic_ztodt, ptend, pbuf)

          ! Destroy massless droplets - Note this routine returns with no change unless
          ! micro_do_massless_droplet_destroyer has been set to true
          call massless_droplet_destroyer(cld_macmic_ztodt, state, ptend)

          ! Limit the value of hydrometeor concentrations in order to place
          ! reasonable limits on hydrometeor drop size and keep them from
          ! becoming too large.
          ! Note:  this needs to be called after hydrometeor mixing ratio
          !        tendencies are adjusted by subcol_SILHS_fill_holes_conserv
          !        and after massless drop concentrations are removed by the
          !        subcol_SILHS_massless_droplet_destroyer, but before the
          !        call to physics_ptend_scale.
          if (trim(subcol_scheme) == 'SILHS') &
            call subcol_SILHS_hydromet_conc_tend_lim( state, cld_macmic_ztodt, ptend )

          ! Copy ptend_aero field to one dimensioned by sub-columns before summing with ptend
          call subcol_ptend_copy(ptend_aero, state_sc, ptend_aero_sc)
          call physics_ptend_sum(ptend_aero_sc, ptend_sc, state_sc%ncol)
          call physics_ptend_dealloc(ptend_aero_sc)

          ! Have to scale and apply for full timestep to get tend right
          ! (see above note for macrophysics).
          call physics_ptend_scale(ptend_sc, 1._r8/cld_macmic_num_steps, ncol)

          call physics_update (state_sc, ptend_sc, ztodt, tend_sc)

          call check_energy_chng(state_sc, tend_sc, 'microp_tend_subcol', &
            nstep, ztodt, zero_sc, &
            prec_str_sc(:state_sc%ncol)/cld_macmic_num_steps, &
            snow_str_sc(:state_sc%ncol)/cld_macmic_num_steps, zero_sc)

          call physics_state_dealloc(state_sc)
          call physics_tend_dealloc(tend_sc)
          call physics_ptend_dealloc(ptend_sc)
        else
          call microp_driver_tend(state, ptend, cld_macmic_ztodt, pbuf)
        end if
        ! combine aero and micro tendencies for the grid
        call physics_ptend_sum(ptend_aero, ptend, ncol)
        call physics_ptend_dealloc(ptend_aero)

        ! Have to scale and apply for full timestep to get tend right
        ! (see above note for macrophysics).
        call physics_ptend_scale(ptend, 1.0_r8 / cld_macmic_num_steps, ncol)

        call diag_clip_tend_writeout(state, ptend, ncol, lchnk, ixcldliq, ixcldice, ixq, ztodt, rtdt)

        call physics_update (state, ptend, ztodt, tend)

        call check_energy_chng(state, tend, 'microp_tend', nstep, ztodt, &
          zero, prec_str(:ncol)/cld_macmic_num_steps, &
          snow_str(:ncol)/cld_macmic_num_steps, zero)

        prec_sed_macmic(:ncol) = prec_sed_macmic(:ncol) + prec_sed(:ncol)
        snow_sed_macmic(:ncol) = snow_sed_macmic(:ncol) + snow_sed(:ncol)
        prec_pcw_macmic(:ncol) = prec_pcw_macmic(:ncol) + prec_pcw(:ncol)
        snow_pcw_macmic(:ncol) = snow_pcw_macmic(:ncol) + snow_pcw(:ncol)
      end do ! end substepping over macrophysics/microphysics

      prec_sed(:ncol) = prec_sed_macmic(:ncol) / cld_macmic_num_steps
      snow_sed(:ncol) = snow_sed_macmic(:ncol) / cld_macmic_num_steps
      prec_pcw(:ncol) = prec_pcw_macmic(:ncol) / cld_macmic_num_steps
      snow_pcw(:ncol) = snow_pcw_macmic(:ncol) / cld_macmic_num_steps
      prec_str(:ncol) = prec_pcw(:ncol) + prec_sed(:ncol)
      snow_str(:ncol) = snow_pcw(:ncol) + snow_sed(:ncol)

    end if

    ! Add the precipitation from CARMA to the precipitation from stratiform.
    if (carma_do_cldice .or. carma_do_cldliq) then
      prec_sed(:ncol) = prec_sed(:ncol) + prec_sed_carma(:ncol)
      snow_sed(:ncol) = snow_sed(:ncol) + snow_sed_carma(:ncol)
    end if

    if (.not. deep_scheme_does_scav_trans()) then

      ! -------------------------------------------------------------------------------
      ! 1. Wet Scavenging of Aerosols by Convective and Stratiform Precipitation.
      ! 2. Convective Transport of Non-Water Aerosol Species.
      !
      !  . Aerosol wet chemistry determines scavenging fractions, and transformations
      !  . Then do convective transport of all trace species except qv,ql,qi.
      !  . We needed to do the scavenging first to determine the interstitial fraction.
      !  . When UNICON is used as unified convection, we should still perform
      !    wet scavenging but not 'convect_deep_tend2'.
      ! -------------------------------------------------------------------------------

      if (clim_modal_aero .and. .not. prog_modal_aero) then
        call modal_aero_calcsize_diag(state, pbuf)
        call modal_aero_wateruptake_dr(state, pbuf)
      endif

      call aero_model_wetdep( state, ztodt, dlf, cam_out, ptend, pbuf)
      call physics_update(state, ptend, ztodt, tend)

      if (carma_do_wetdep) then
        ! CARMA wet deposition
        !
        ! NOTE: It needs to follow aero_model_wetdep, so that cam_out%xxxwetxxx
        ! fields have already been set for CAM aerosols and cam_out can be added
        ! to for CARMA aerosols.
        call carma_wetdep_tend(state, ptend, ztodt, pbuf, dlf, cam_out)
        call physics_update(state, ptend, ztodt, tend)
      end if

      call convect_deep_tend_2(state,   ptend,  ztodt,  pbuf)
      call physics_update(state, ptend, ztodt, tend)

      ! check tracer integrals
      call check_tracers_chng(state, tracerint, 'cmfmca', nstep, ztodt,  zero_tracers)
    end if

    !===================================================
    ! Moist physical parameteriztions complete:
    ! send dynamical variables, and derived variables to history file
    !===================================================

    call diag_phys_writeout(state, pbuf)
    call diag_conv(state, ztodt, pbuf)

    !===================================================
    ! Write cloud diagnostics on history file
    !===================================================

    call cloud_diagnostics_calc(state, pbuf)

    !===================================================
    ! Radiation computations
    !===================================================

    call radiation_tend(state, ptend, pbuf, cam_out, cam_in, net_flx)

    ! Set net flux used by spectral dycores
    do i = 1, ncol
      tend%flx_net(i) = net_flx(i)
    end do

    call physics_update(state, ptend, ztodt, tend)

    call check_energy_chng(state, tend, 'radheat', nstep, ztodt, zero, zero, zero, net_flx)

    ! Diagnose the location of the tropopause and its location to the history file(s).
    call tropopause_output(state)

    ! Save atmospheric fields to force surface models
    call cam_export(state, cam_out, pbuf)

    ! Write export state to history file
    call diag_export(cam_out)

  end subroutine tphysbc

  subroutine phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d)

    !-----------------------------------------------------------------------------------
    !
    ! Purpose: The place for parameterizations to call per timestep initializations.
    !          Generally this is used to update time interpolated fields from boundary
    !          datasets.
    !
    !-----------------------------------------------------------------------------------

    use shr_kind_mod        , only: r8 => shr_kind_r8
    use chemistry           , only: chem_timestep_init
    use chem_surfvals       , only: chem_surfvals_set
    use physics_types       , only: physics_state
    use physics_buffer      , only: physics_buffer_desc
    use carma_intr          , only: carma_timestep_init
    use ghg_data            , only: ghg_data_timestep_init
    use cam3_aero_data      , only: cam3_aero_data_on, cam3_aero_data_timestep_init
    use cam3_ozone_data     , only: cam3_ozone_data_on, cam3_ozone_data_timestep_init
    use aoa_tracers         , only: aoa_tracers_timestep_init
    use vertical_diffusion  , only: vertical_diffusion_ts_init
    use radheat             , only: radheat_timestep_init
    use solar_data          , only: solar_data_advance
    use qbo                 , only: qbo_timestep_init
    use iondrag             , only: do_waccm_ions, iondrag_timestep_init
    use prescribed_ozone    , only: prescribed_ozone_adv
    use prescribed_ghg      , only: prescribed_ghg_adv
    use prescribed_aero     , only: prescribed_aero_adv
    use aerodep_flx         , only: aerodep_flx_adv
    use aircraft_emit       , only: aircraft_emit_adv
    use prescribed_volcaero , only: prescribed_volcaero_adv
    use prescribed_strataero, only: prescribed_strataero_adv
    use mo_apex             , only: mo_apex_init
    use epp_ionization      , only: epp_ionization_active

    type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
    type(cam_in_t),      intent(inout), dimension(begchunk:endchunk) :: cam_in
    type(cam_out_t),     intent(inout), dimension(begchunk:endchunk) :: cam_out

    type(physics_buffer_desc), pointer :: pbuf2d(:,:)

    ! Update geomagnetic coordinates
    if (epp_ionization_active .or. do_waccm_ions) then
      call mo_apex_init(phys_state)
    end if

    ! Chemistry surface values
    call chem_surfvals_set()

    ! Solar irradiance
    call solar_data_advance()

    ! Time interpolate for chemistry.
    call chem_timestep_init(phys_state, pbuf2d)

    ! Prescribed tracers
    call prescribed_ozone_adv    (phys_state, pbuf2d)
    call prescribed_ghg_adv      (phys_state, pbuf2d)
    call prescribed_aero_adv     (phys_state, pbuf2d)
    call aircraft_emit_adv       (phys_state, pbuf2d)
    call prescribed_volcaero_adv (phys_state, pbuf2d)
    call prescribed_strataero_adv(phys_state, pbuf2d)

    ! Prescribed aerosol deposition fluxes
    call aerodep_flx_adv(phys_state, pbuf2d, cam_out)

    ! CAM3 prescribed aerosol masses
    if (cam3_aero_data_on) call cam3_aero_data_timestep_init(pbuf2d,  phys_state)

    ! CAM3 prescribed ozone data
    if (cam3_ozone_data_on) call cam3_ozone_data_timestep_init(pbuf2d,  phys_state)

    ! Time interpolate data models of gasses in pbuf2d
    call ghg_data_timestep_init(pbuf2d, phys_state)

    ! Upper atmosphere radiative processes
    call radheat_timestep_init(phys_state, pbuf2d)

    ! Time interpolate for vertical diffusion upper boundary condition
    call vertical_diffusion_ts_init(pbuf2d, phys_state)

    !----------------------------------------------------------------------
    ! update QBO data for this time step
    !----------------------------------------------------------------------
    call qbo_timestep_init()

    call iondrag_timestep_init()

    call carma_timestep_init()

    ! age of air tracers
    call aoa_tracers_timestep_init(phys_state)

  end subroutine phys_timestep_init

end module physpkg
